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Abstract: We study the impact of initial-state multiple parton radiation on transverse 
momentum (qx) distribution of Higgs boson produced via bottom quark fusion at hadron 
colliders. The shape of the resulting qx distribution is affected by the bottom-quark mass 
corrections and by the strong kinematical behavior of the bottom-quark parton density. We 
account for both features in the full range of qr- To do this, we formulate the resummation 
calculation in a general-mass factorization (S-ACOT) scheme and introduce a correction in 
the resummed-term to account for the effect from large-gr kinematics of Higgs boson. The 
results of this resummation are compared to fixed-order and PYTHIA predictions. 
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1. Introduction 

Understanding of electroweak symmetry breaking (EWSB) is the central challenge for high 
energy physics. The search for Higgs boson(s), assumed to be responsible for the generation 
of gauge-boson and fermion masses, is the primary task for the existing and future high 
energy colliders. 

The Higgs sector may be represented by one complex scalar doublet, as it is economi- 
cally realized in the Standard Model (SM), or by two or more doublets, as it takes place in 
the Minimal Supersymmetric Standard Model (MSSM) and its extensions. In MSSM, two 
Higgs-doublet superfields are necessary to generate masses for up- and down-type quarks 
and provide remarkable cancellation of triangle anomalies. 

The two MSSM Higgs doublets have independent vacuum expectation values (VEVs), 
v u and Vd- The sum of the squares of these VEVs is fixed by the well-known Z boson mass, 
but their ratio, denoted as tan/3 = v u /vd, is a free parameter of the model. As a result of 
spontaneous symmetry breaking, five physical particles appear in the Higgs sector: h (light 
CP-even), H (heavy CP-even), A (CP-odd), and if ± (charged). An important feature of 
MSSM is that, for large values of tan/3, the Yukawa couplings of the 6-quarks to the neutral 
Higgs boson TL (where TL = h, H, or A) are enhanced by a factor 1/ cos/3 compared to the 
SM bbHsM Yukawa coupling. 

In the MSSM, Higgs boson masses are functions of CP-odd Higgs mass, tua, and tan/3. 
The experimental lower limit on the Higgs boson mass deduced from the LEP data [1] 
favors scenarios with a intermediate-to-large values of tan/3 (> 5 — 10). Theoretically, 
scenarios with large tan/3 are highly motivated by SO(10) models of supersymmetric grand 
unification [2-19]. Therefore, processes involving Higgs boson production via enhanced bbTL 
Yukawa coupling in the MSSM and other new physics models may have large production 
rates and could play an important role in the study of the Higgs boson [20]. 

The partonic processes contributing to the inclusive Higgs boson production with en- 
hanced bbTL coupling are represented by (a) bb — ► TL; (b) gb — ► TLb; and (c) gg — ► bbTL 
scattering. In perturbative quantum chromodynamics (QCD), these are not independent 
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production mechanisms, since b partons inside the hadron beam/target arise from QCD 
evolution (splitting) of gluons, and gluons radiate off quarks [21-23]. The three processes 
(a,b,c) all give rise to the same hadronic final states, with two S-mesons appearing in dif- 
ferent, but overlapping, regions of phase space — either as beam/target remnants, or as high 
transverse momentum particles. The distinction between the three processes depends very 
much on the factorization scheme adopted for the QCD calculation, as has been recently 
reviewed in Ref. [20]. 

These TLbb processes have been extensively studied recently in SM and MSSM scenar- 
ios [24-37]. Calculations in the 5-fiavor scheme have been carried out to the 2-loop level 
[32], which considerably reduced the theoretical uncertainty due to the perturbative expan- 
sion, as estimated by the residual scale dependence. Comparison of results obtained in the 
4- and 5-fiavor schemes has also been carried out [36, 37]. It shows consistency between the 
two schemes in the energy region of the Fermilab Tevatron and the CERN Large Hadron 
Collider (LHC). 

In spite of the good theoretical control of the inclusive Higgs boson production cross 
section in bb fusion, neither the 4-fiavor nor 5-fiavor scheme is adequate for predicting the 
transverse momentum (qx) distribution of the Higgs boson, when qT values are of the order 
of the b quark mass (rrib), at any fixed order of the perturbative calculation. To properly 
describe the low-gr region and determine the qx range where fixed-order calculations are 
applicable, one should resum large soft and collinear logarithms, as we discuss in detail 
in the next section. Furthermore, a more comprehensive factorization scheme (general- 
mass variable flavor number scheme) must be used to describe non-negligible dependence 
on m;, when qx is comparable to m^. Many studies of soft gluon resummation for the 
kinematical distributions of Higgs boson produced via gluon-gluon fusion are available in 
the literature [38-45]. As for the soft gluon resummation in production of Higgs boson 
via bb fusion, it was first discussed for massless b quarks in Ref. [46]. Later, Ref. [35] has 
studied the effect of 0(a s ) subleading logarithmic contribution in the Sudakov form factor, 
but only in the soft-gluon limit. 

As shown in Refs. [47,48], the correct model for the transverse momentum distribution 
of the Higgs boson is crucial for unambiguous reconstruction of Higgs boson mass in the 
TL — > tt decay channel. It is also important for discriminating the signal events from the 
backgrounds by carefully examining the qT distribution of the Higgs boson in TLbb associ- 
ated production, followed by TL — ► bb decay [49]. In this work, we study the effect of the 
initial-state multiple parton radiation on the transverse momentum distribution of Higgs 
boson produced at hadron colliders via bb fusion. Our calculation includes contributions 
from resummation up to next-to-next-to-leading logarithmic accuracy and fixed-order per- 
turbation theory up to next-to-leading order. 1 In view of the non-negligible mass of the 
bottom quark and strong kinematical dependence of the 6-quark parton density, special 
corrections (referred in this work as the "heavy-quark mass correction" and "kinematical 
correction") must be included in the resummation formalism to predict the bb — > TL cross 
section in the full range of qx- 

lr This order of accuracy is the same as that in the calculation [50] of the kinematical distributions for 
weak gauge bosons produced via light quark annihilation in hadron collisions. 
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The rest of this paper is organized as follows. In Section 2, we discuss the theoretical 
framework for transverse momentum resummation used in this study and introduce the 
"heavy-quark mass correction" and "kinematical correction'". In Section 3, we present the 
numerical results for bb — > H at the Tevatron and the LHC and compare the qr distributions 
predicted by resummation calculations to PYTHIA predictions. Section 4 contains our 
conclusions. 

2. Resummation for heavy flavors: theory 

Kinematical distributions of Higgs boson produced in bottom quark-antiquark fusion at 
hadron colliders are often computed within the zero-mass variable flavor number (ZM-VFN) 
factorization scheme, which neglects masses of bottom and lighter quarks in hard-scattering 
amplitudes. The 5-flavor scheme is a realization of the ZM-VFN scheme at energies above 
the bottom-quark mass threshold. The application of this scheme is justified because the 
mass M-h of the Higgs boson of interest lies in the range of a few hundred GeV, i. e. , it is much 
larger than the mass nib of the bottom quark. If all momentum scales in the cross section are 
of order Q ~ 3> m^, for instance, when the inclusive rate for Higgs boson production is 
computed, the mass dependence of the hard-scattering amplitude is suppressed by powers of 
a small factor m 2 /Q 2 and can be safely discarded. Meanwhile, large logarithms ln p (Q 2 /m 2 ), 
arising from bb pairs being produced in collinear splittings of gluons, are resummed in the 
parton distribution function (PDF) of 6-quark in the proton by utilizing Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) equations [51-54]. The perturbative parton distribution 
of the bottom quark is set identically equal to zero at factorization scales fxp below the 
6-quark mass threshold (/j,f < fUb) and turned on at /j,f > rub- In technical terms, such 
a definition corresponds to regularization of the ultraviolet divergence in the perturbative 
6-quark PDF by zero-momentum subtraction at hf < rrib and MS subtraction at fj,p > "m-b 
[55]. 

The 4-flavor scheme can be viewed as a realization of the fixed-flavor number, or FFN, 
scheme [56-61], as far as the bottom quarks are concerned. This scheme does not introduce 
nontrivial parton densities for heavy quarks, but rather keeps the heavy-quark contributions, 
including the logarithmic terms \n p {Q 2 / m 2 ) , in the hard-scattering amplitudes. The heavy- 
quark PDF is set identically to zero in the FFN scheme, as a consequence of regularization 
by zero-momentum subtraction in the whole range of /if- As in the ZM-VFN scheme, 
the mass-dependent terms are negligible in the FFN scheme at Q 2 3> m 2 , except for the 
logarithms ln p (Q /m^), which, however, are not resummed into the heavy-quark parton 
densities. Hence, at large Q, the factorization scale dependence in the FFN perturbative 
calculation is stronger than that in the ZM-VFN calculation, though the FFN calculation 
is more reliable when Q is close to mb. Both ZM-VFN and FFN schemes were applied 
recently to Higgs boson production via bb fusion with zero, one or two tagged b quarks in 
the final state [27,28,31-33,36,37]. A recent computation of 0(a 2 ) corrections [32] to the 
bb — > 7i cross section was realized in the ZM-VFN scheme, since implementation of the 
quark mass dependence is tedious beyond 0(a s ). 
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However, neither the ZM-VFN scheme nor the FFN scheme will work well in a calcu- 
lation of the transverse momentum distribution of the Higgs boson, when qr is small, of 
the order of the bottom quark mass. To properly describe the small-gr region, large soft 
and collinear logarithms ln m (q T /Q 2 ) must be resummed together with the mass logarithms 
ln p (Q 2 /m^), while keeping the essential dependence on the 6-quark mass at qx ~ "mb- 
Soft gluon resummation must be realized in a more comprehensive factorization scheme (a 
general-mass variable flavor number scheme [62]) to achieve this goal [63]. Formulation of 
the Collins-Soper-Sterman (CSS) resummation method [64] in such a scheme (simplified 
Aivazis-Collins-Olness-Tung, or S-ACOT, scheme [62,65]) has been proposed recently [63] 
to compute qx distributions in scattering processes initiated by heavy quarks. The ad- 
vantage of the S-ACOT scheme is that it retains massless expressions for hard-scattering 
amplitudes with incoming heavy quarks (flavor-excitation amplitudes), and hence, drasti- 
cally simplifies the calculation, while preserving the relevant mass terms. In this work, we 
adopt the qr resummation formalism with the heavy-quark (HQ) mass effects proposed in 
Ref. [63], hereafter referred to as "the CSS-HQ formalism". 

The CSS formalism and its application to bb — > TC in the ZM-VFN scheme has been 
discussed extensively in Refs. [46,50,66], and we refer the reader to those publications for 
details. Here we summarize the features of the CSS formalism essential for the understand- 
ing of the numerical results in the next section. We symbolically write the total (TOT) 

resummed differential cross section da/dQ 2 dydq T for h\(P\)h2(P2) T~t(q)X, where h± 
and \i2 are the initial-state hadrons, as 



W denotes the Fourier-Bessel integral of a resummed form factor W(b, Q, x\,X2) introduced 
in space of the impact parameter b (Fourier-conjugate to qx): 



PERT is the perturbative QCD cross section, evaluated at a finite order of the QCD coupling 
strength a s . The asymptotic piece ASY is defined as a perturbative QCD expansion of the 
resummed W-term to the same order of a s as in PERT. The difference of PERT and ASY, 
which we denote by Y = PERT— ASY, is the regular part of the cross section. x\ and X2 are 
the light-cone momentum fractions of the partons entering the WH vertex in W and ASY. 
They satisfy x lj2 = Qe ±v j\[S at q T -> 0, with Q and y = \ ln((g° + q 3 )/{q° - q 3 )) being the 
invariant mass and rapidity of the Higgs boson, and the square of the h\h2 center-of-mass 
energy S = (Pi + P2) 2 ■ We will discuss the choice of xi j2 at nonzero qT in a moment. 

The resummed form factor W(b, Q, x±,x 2 ) is composed of perturbative QCD contribu- 
tions (dominant at b — ► 0) and non-perturbative contributions (dominant at b — > 00). We 
adopt the b* prescription [64,67] to evaluate the form factor W(b, Q, x±, X2) numerically, 
which is given by 



TOT = W + PERT - ASY. 



(2.1) 




W(b,Q, Xl ,x 2 ) = W pert (b*,Q,x 1 ,X2)e- :F » p ( b > Q \ 



(2.2) 
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In Eq. (|2.2D , W per t(b*, Q,xi,x 2 ) is the perturbative form factor, evaluated as a function of 



the variable 6* = bj a/1 + b 2 /6^ a2 , with 6 max being a free parameter: 

ai,02 

[(Cfc/ 0l ® fa 1 /h 1 )(xi,b) (C- b/a2 f a2/h2 )(x 2 ,b) + (b 



b) 



(2.3) 
(2.4) 



Here 



S P (b,Q) 



c ^ 2 dtf 



A(a,{fi))hx[ ^ ) + £(a s (/2)) 



(2.5) 



is the perturbative Sudakov factor, with bo = 2e~ lE ~ 1.123, and 

rl x 
^2(Cb/a® fa/h)(x,b) = -rC b/a ^,b,fi F ,m b )f a/h (-,fj, F ) (2.6) 

is a convolution of the Wilson coefficient function C b i a (x,b, //p,m&) for outgoing b quarks 
and parton density f a /h( x i ^f), summed over the intermediate parton states a. The overall 
normalization is given by 



r (o) 



7r f m b (Q) tan/3 
6 



5(Q 2 - M%) , 



(2.7) 



where f3 and v are the Higgs-doublet mixing angle and vacuum expectation value, respec- 
tively. The bottom quark running mass mb(Q) is evaluated at the scale Q [46], and is 
the mass of the Higgs boson H. 

The non-perturbative contributions are described by the exponential e ~~ FNP ^ b ' < ^\ For 
simplicity, the function J r /vp(6, Q) is assumed to be the same in the heavy- and light-quark 
channels and taken from the global q F fit [66] for b max = 0.5 GeV -1 . This approximation is 
sufficient for the purposes of the comparison of the resummation calculation with PYTHIA 
and may be refined by using alternative prescriptions, such as the revised 6* model [68], in 
the future analyses. 

The perturbative expressions for the functions PERT, ASY, and the coefficients A- 1 ', 
an d ,4(2) of the f unct i n W can be found in Refs. [46,50]. The 0(a s /ir) coefficients 

C^)(x, b, /J,p) in the Wilson coefficient functions C b / a (x, b, /j, f ) are different in the ZM-VFN 

and S-ACOT factorization schemes. In the ZM-VFN scheme, 



cff b (x,b,fJ, F ) = ~y(1-x) - 
+ C F 5{l-x) 
& ' V>f) = T R x(l - x) 



-In 
In 



°±e-»'A + - + — 
b C 2 ) 4 16 



(2.8) 
(2.9) 



- 5 - 



where T R = 1/2 and C F = (N% - l)/(2iV c ) = 4/3 are the QCD Casimir invariants (with 
N c = 3), Pff q {x) = C F [(1 + x 2 )/(l - x)}+ and = ^(x 2 + (1 - x) 2 ) are the 0{a s ) 

splitting functions, and V =ir 2 — 2. 2 

The coefficient for Nf active quark flavors can be obtained following the method 
presented in Refs. [69,70]. It is found to be 



with fa = (HJVc - 2JV/)/6, 



universal + A) J (2.10) 



8£L- = + ^ + AO, ( (in £) 2 - I In C 2 - ^ ft) 

-*((S-3)*-S^)-(f)- 



and 



^.-^f|-T + «&V^('S + ^-36 , )-^/('i + ^'l ■ (2-12) 



where C3 = 1.202057... is the Riemann zeta-function C, n for n = 3. Throughout the paper, we 
choose the factorization constants C\ = bo and C2 = 1, and the factorization scale /xp = 60/6 
in the Wilson coefficient functions C b / a (x,b, fj,p) and parton distributions f a /h{ x i A*f)- We 
note that the above results for the and £>( 2 ) coefficients are different from the ones 
presented in Ref. [35], where only soft gluon contributions were considered. 

In the S-ACOT scheme, we include m b dependence in the Wilson coefficient C^)(x, b, fj,p) 
for gluon splitting into a bb pair (i.e., b <— g splittings): 

cly g (x,b,m b ,^ F ) = T R x(l - x)bm b Ki(bm b ) 

M\ / II T7>\ 

(2.13) 



K (bm b ) - 9(hf ~ m b )ln(—) 

\m b J 



where Kq(z) and Ki(z) are the modified Bessel functions [71]. This expression reduces to 
the massless result in Eq. fl2.9| ) when b <C m b : 

lim 

bnrih— 



C§ g (x,b,m b ,VF) = T R x(l - x) - In f^j P^(x) , (2.14) 



since Kq(z) — ► and K\(z) — > I/2 as z -* 0. We keep massless expressions for the 
remaining perturbative coefficients in W, PERT, and ASY, in accordance with the rules of 
the S-ACOT scheme [63]. In particular, the mass dependence is neglected in PERT and 
ASY, as their difference contributes substantially to Eq. ( p] ) only at qx ~ Q m b . 

While PERT is evaluated by using the exact 2^2 kinematics, the phase space element 
in the W and ASY terms is unique only in the limit qx — ► 0, but may be defined in 



2 The V term shown here is different from that in Eq. (14) in Ref. [46], as a result of a typo in Eq. (14), 
cf. Eq. (A2), of that paper. 
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several ways at nonzero qr- At qr of order Q, momentum fractions x of the initial-state 
partons must be large enough to produce the Higgs boson with a large transverse mass 
Mt = \J Q 2 + Or'' too sma ll momentum fractions x ~ Q/\/~S accessible at qr —> are not 
allowed. To introduce information about the reduction of phase space available for collinear 
QCD radiation at large qr, we define the light-cone momentum fractions x\ and X2 in W 
and ASY as 

xi, 2 = 2(P 2) l • q)/S = M t e ±y /VS . (2.15) 

As qr — > 0, x\ and x<i reduce to their Born-level expressions, xP° rn = Qe ±y /\^S, and 
the canonical CSS form is reproduced. At larger qr, contributions from unphysical small 
momentum fractions are excluded from W and ASY by the growing M t in x\p. We will 
argue that this prescription of redefining the values of x\ t 2, to be referred as the "kinematical 
correction" in this work, is crucial for predicting the resummation cross sections in bb channel 
at large qr, due to the strong dependence of the 6-quark PDF on x in the relevant kinematic 
regions of the Tevatron and the LHC. 

The next section presents numerical comparisons of the parts of the resummation cross 
sections (TOT, W, PERT, and ASY), computed under various assumptions about the order 
of the perturbative coefficients in a s , factorization scheme, and kinematical correction. Our 
naming conventions for various terms are summarized in Table |IJ A numerical argument in 
the names of the PERT and ASY functions indicates the order of the perturbative calcu- 
lation. For example, PERT(l) stands for the 0(a s ) perturbative piece. Three arguments 
following the names of the functions W and TOT indicate orders of the functions A and B 
in the Sudakov factor ( |2,5| ), and Wilson coefficient functions C b / a (x,b, fiF,i^b) in Eq. Q2.6| ). 
For example, the .A, £>, and C functions in the W(2,2,l) distribution are evaluated up to or- 
ders a 2 s , a 2 ., and a s , respectively. The S-ACOT factorization scheme and active kinematical 

correction are indicated by the superscript "HQ" (for "heavy-quark") and subscript "KC" (for 

HO 

"kinematical correction"), as in W^(1,1,0) for the W(l, 1, 0)-term with the kinematical 
correction evaluated in the S-ACOT scheme. 





Order of QCD coupling strength a s 


Kinematical 


S-ACOT 




A{a s (p)) 


B(a,(fi)) 


Cb/a 


Y 


correction 


scheme 


TOT(1,1,0) 


1 


1 





1 






TOT(2,2,l) 


2 


2 


1 


1 






TOT KC (2,2,l) 


2 


2 


1 


1 


yes 




W(1,1,0) 


1 


1 











W(2,2,l) 


2 


2 


1 








W K c(2,2,l) 


2 


2 


1 




yes 




W^(2,2,l) 


2 


2 


1 




yes 


yes 


PERT(l) 








1 






ASY(l) 








1 






ASY KC (1) 








1 


yes 





Table 1 : Naming conventions for qr distributions in Section ||. 



-7- 



3. Numerical results 



In this section we present the numerical results of our study for a Higgs boson produced 
via bottom quark fusion at the Tevatron Run-2 (a 1.96 TeV proton-antiproton collider) and 
the LHC (a 14 TeV proton-proton collider). To make our study less dependent on SUSY 
parameters, we focus on production of the CP-odd Higgs particle A, since the Yukawa 
couplings of A to the heavy quarks are independent of the Higgs mixing angle. 

First, we study in details the Tevatron case, and then present our final results also for 
the LHC. As an example, we have chosen uia to be 100 GeV and 300 GeV, respectively, for 
the Tevatron and LHC studies. Our results could be trivially generalized to the production 
of any neutral Higgs boson 7i by properly scaling the bbA Yukawa coupling to the actual 
one. In our studies the value of tan/3 was chosen to be 50. The tree level Yukawa coupling of 
bottom quarks and the CP-odd Higgs boson (A) in the minimal supersymmetric Standard 
Model (MSSM) is equal to m& tan (3/{y2v) , where the bottom quark mass is 4.7 GeV, and 
the vacuum expectation value v is about 246 GeV. To improve the numerical predictions, we 
have resummed the large logarithms originated from the QCD radiative corrections to the 
bbA Yukawa coupling by introducing the running bottom quark mass at the scale of Higgs 
boson mass, as done in Ref. [46] for calculating the next-to-leading order QCD corrections 
to the production of bb — > A in hadron collisions. The running bottom quark mass is about 
2.98 GeV for m A = 100 GeV. 

The full event generator, such as PYTHIA, can give a fair description of the event 
topology after turning on the QCD showering from the initial and, possibly, final states, 
which are generated by the probability functions calculated from the relevant Sudakov form 
factors [72]. Therefore, it is a common practice to compare PYTHIA predictions to the 
event shape of experimental data. For example, the shape of the qx distribution of the vector 
boson produced at hadron colliders can be fairly described by the PYTHIA program, though 
the normalization of the event rate is usually off the scale, because it does not include all 
the finite part of higher order QCD corrections [50]. As compared to the CSS resummation 
formalism, the PYTHIA calculation does not include contributions generated from the C- 
functions and the Y-term. Therefore, a fair comparison to the PYTHIA prediction is to 
include only the Sudakov contributions in a resummation calculation. In Fig. |l|, we show the 
prediction of W(1,1,0) and compare it to the PYTHIA prediction for the hard scattering 
process bb — > H [72]. 3 It is evident that W(1,1,0) predicts a very different shape of qx 
distribution from PYTHIA in bb — > A production, though the integrated rates (i.e., the 
areas under the two curves) are about the same. On the contrary, in the case of vector boson 
production that is dominated by light quark scatterings in the initial state, the above two 
calculations predict similar, though not identical, qx distributions [50]. Fig. |] illustrates 
that the peak position of qx distribution predicted by PYTHIA is lower than that by the 
resummation calculation W(1,1,0). Also, PYTHIA predicts a narrower shape in the qx 
distribution. It is important to understand the cause of this difference in order to reliably 

3 The invariant mass cutoff of parton showers (the parameter PARJ(82) in the computer code), below 
which partons are assumed not to radiate, was chosen to be 1.0 GeV, which is the default value of PYTHIA 
version 6.2. 
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predict the qx distribution of the Higgs boson produced in bb — > A process. 

A close examination reveals that when PYTHIA generates QCD showering, the kine- 
matical distributions of the final-state particles (including the quarks and gluons generated 
from QCD showering) are slightly modified to satisfy energy-momentum conservation at 
each stage of showering. Namely, the kinematics of Higgs boson is modified according to 
the amount of showering. By this, it effectively includes some part of higher-order con- 
tribution (similar to part of PERT(l) contribution). On the other hand, in the W(1,1,0) 
calculation the emitted soft gluons are assumed not to carry any momentum (strictly in the 
soft limit), which implies that the Y-term contribution could be important for determining 
the shape of qx distribution. From Fig. [j], we expect the Y-term contribution to be nega- 
tive at large qx ~ rriA in order for the result of resummation calculation to resemble more 
the PYTHIA prediction. Moreover, the Sudakov form factors in PYTHIA contain some 
additional 0{a 2 s ) contributions when using the next-to-leading order PDF's to generate 
the event distributions. Thus, we should improve the resummation calculation by includ- 
ing next-to-leading order perturbative corrections. The resulting distribution, denoted as 
TOT(2,2,l), is the sum of W(2,2,l) and PERT(l) with the subtraction of ASY(l) to avoid 
the overlapped contribution, cf. Eq. (fO|). 

As shown in Fig. |, ASY(l) and PERT(l) almost coincide as q T -> 0, and then ASY(l) 
curve goes increasingly above PERT(l) curve as qT increases, yielding a negative term 
Y=PERT(1)-ASY(1). By adding such a negative Y-term to W(2,2,l) in TOT(2,2,l), we 




Figure 1: qx distributions predicted by W(l, 1, 0) and PYTHIA, shown as solid and dashed curves, 
respectively, for a 100 GeV Higgs boson produced via bb fusion at the Tevatron Run-2. 
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improve agreement of the resummation and PYTHIA predictions. This behavior confirms 
our expectations from the above mentioned inspection of the distributions shown in Figjl]. 
It is also important to know what causes the difference from the vector boson produc- 
tion in light-quark annihilation channels (e.g., uu — > Z), in which the Y-term is positive. 
Two possible causes are the differences between the hard-part matrix elements included 
in PERT(l) for bb — > A and uu — > Z processes, and the differences in the shapes of the 
parton distributions for b and u quarks at the x values typical for production of a 100 GeV 
Higgs boson at the Tevatron. This will be illustrated in the next figure. Another comment 
about Fig. I is that the difference in TOT(2,2,l) and TOT(1,1,0) comes from the difference 
in W(2,2,l) and W(1,1,0). The coefficients included in W(2,2,l) mainly change the 
overall normalization to include (part of) the next-to-leading order contributions that are 
not accounted for, even after including higher-order (i.e., 0(a 2 s )) contributions and 
£>( 2 ) in the Sudakov form factor. 

To examine the effects of the hard-part matrix elements and the shape of parton dis- 
tribution functions on the shape of qr distribution, Fig. || compares the results of four 
calculations. The upper-left plot shows again the result of our bb — > A calculation pre- 
sented in Fig. ||. The upper-right plot is obtained from the upper-left plot by replacing 
the hard-part matrix element PERT(l) for producing a CP-odd Higgs boson A by that for 
producing a vector boson Z. (Note that the quark masses are neglected in both cases.) The 




Figure 2: qr distributions contributed by various pieces of the resummation calculations. ASY(l) 
and PERT(l) almost coincide as qr — > 0, and ASY(l) curve goes increasingly above PERT(l) curve 
as <7t increases to yield a negative term Y=PERT(1)-ASY(1). 
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overall normalization factor contributed by the bbA vertex is not changed. Consequently, 
the difference between the two upper plots is entirely due to the kinematical dependence in 
the hard-part matrix elements in the two processes. 

The lower-left plot is obtained from the upper-left plot for bb — > A by replacing the b 
and b PDF's by u and u PDF's. This is to examine the effect from the different shapes 
of the b- and u-quark PDF's. Finally, the lower-right plot is obtained from the upper-left 
plot by replacing both the hard-part matrix element for bb — > A and 6-quark PDF by the 
hard-part matrix element for uu — > Z and w-quark PDF. 

The four plots clearly show that the difference due to the replacement of the hard-part 
matrix element is small, while the difference due to the replacement of the PDF is large. 
If we replace the bottom quark PDF by the up quark PDF, then the Y-term becomes 
positive, in agreement with what has been observed in the previous calculations for vector 
boson production. From this comparison, we conclude that the shape of bottom quark PDF 
causes the Y-term to be negative. As all sea-parton PDF's, the 6-quark PDF is a rapidly 
decreasing function of x in the whole range of x, while the u-quark PDF, which includes a 




Figure 3: Sensitivity of qr distributions to the choice of hard-part matrix elements in PERT(l) 
and parton distribution functions. The upper-left plot is the reproduction of Fig. 2. In the upper- 
right plot, the hard-part matrix element of bb ~ * A is replaced by that for uu — > Z, but the PDF 
is not changed. In the lower-left plot, the bottom-quark PDF is replaced by the up-quark PDF, 
but the hard-part matrix element of bb — ► A is not changed. In the lower-right plot, both the Z 
hard-part matrix element and u-quark PDF are used in the calculation. The plots show that the 
Y-term in bb — > A is negative mainly due to the strong x dependence of the bottom quark PDF. 
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substantial valence component, varies slower with x. As discussed in Section 2, the energy- 
momentum conservation in PERT(l) prevents too small x Q/V~S from contributing at 
nonzero qx- Meanwhile, such x are included in ASY(l), evaluated using small-gr phase 
space. Consequently ASY(l) may be enhanced compared to PERT(l) by large PDF's from 
smaller x, especially in the case of the sea-quark PDF's. As a result of the PDF-induced 
enhancement in bb scattering, inclusion of the Y-term brings the resummation calculation 
closer to the PYTHIA prediction. 

As discussed in Section 2, the "kinematical correction" (KC), which compensates for 
small, but nonzero energy of soft gluon emissions, can play an important role in explaining 
the difference between the resummation and PYTHIA predictions. Thus, we calculated 
TOT KC (2,2,l) by adding W KC (2,2,1) and PERT(l) and subtracting ASY KC (1) to compare 
it with TOT(2,2,l) in Fig. ||. It is indeed the case that the kinematical correction reduces 
the rate in the large-^T region, which is consistent with our intuition that qr will be re- 
duced by the emission of (soft) gluons carrying away some energy. Fig. ^| also shows that 
the location of the peak of the q? distribution is not affected by the kinematical correction. 
Hence, to explain why the peak position predicted by PYTHIA is lower than that in the 
resummation calculations discussed above, we should take into account another important 
piece of physics, which is included in the PYTHIA event generator, but not in the above 
resummation calculations. That is the effect of the bottom-quark mass mi on the uninte- 
grated parton distribution of the bottom quark inside the resummed W term. PYTHIA 




Figure 4: Impact of the kinematic correction on the qr distribution. The kinematical correction 
reduces the rate in the large-qr region. 
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always generates the initial-state bottom quarks from the gluon splittings, with the proper 
kinematics for the massive bb pairs imposed at the end of the backward-radiation cascade 
chain [72]. However, the proper mj, dependence was not included in the massless resum- 
mation calculation discussed thus far. The mj dependence can be introduced in the CSS 
resummation formalism using the method developed in Ref. [63] and discussed in detail in 
Section 2. The effect of the heavy-quark mass correction can be clearly seen from Fig. |[. 
While the kinematical correction reduces the rate at large qx, the heavy mass m& shifts the 
peak of the distribution to lower qx, so that the resummation cross section agrees better 
with the PYTHIA cross section. 

In conclusion, our prediction on the qx distribution of Higgs boson produced via bot- 
tom quark fusion in hadron collisions is given by TOT(l), which is obtained by adding 
W^q (2,2,1) and PERT(l) and subtracting ASYkc(I)- The numerical result is shown in 
Fig. | (left) for Tevatron (m A = 100 GeV) and in Fig. | (right) for LHC (m A = 300 GeV), 
respectively, where TOT(l) is also compared to the PYTHIA prediction and the fixed-order 
prediction PERT(l). 

As one can see, results for Tevatron and LHC are qualitatively similar. It remains to 
be the case that the peak position in qx distribution predicted by PYTHIA is lower than 
that by TOT(l), and in the large qx region, TOT(l) rate is larger than the PYTHIA rate. 




5 10 15 20 25 30 35 40 45 50 

qjHiggsfGeV) 



Figure 5: The combined effect of the heavy-quark mass (HQ) and kinematical (KC) corrections. 
The heavy-quark mass correction shifts the peak position of the qr distribution to a lower qx , while 
the kinematical correction reduces the rate in the large-qr region. 
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Next, we examine how our conclusions are modified by the 0(a 2 ) corrections to the Y- 
term. We calculate ASYkc(2) using its exact 0(a 2 ) expression and estimate the PERT(2) 
piece in the large-gy region by multiplying PERT(l) by the 0(a 2 ) "K-factor" (~ 1.75) for 
the Tevatron case, extracted from Fig. 11 in the recent calculation [27] of the 0{a 2 ) rate for 
Higgs boson production at large qx- The resulting TOT(2), evaluated similarly to TOT(l), 
but using the exact ASYkc(2) and the estimated PERT(2), is only slightly smaller than 
TOT(l), not more than a few percent, for qx larger than 30 GeV for the Tevatron case. 

Finally, Fig. [?] shows the integrated cross section as a function of the minimal qx in the 
calculation for the Tevatron (left) and LHC (right). This is another way to illustrate the 
differences in the shape of qr distribution obtained in the resummation, fixed-order, and 
PYTHIA calculations. 

4. Discussion and Conclusion 

We studied the effect of initial-state multiple soft-gluon radiation on the transverse mo- 
mentum (c/t) distributions of Higgs boson produced via bottom-quark fusion at hadron 
colliders. Due to the shape of the bottom-quark parton distribution function that rapidly 
decreases with x, the Y-term in the bb — > A process is negative, and the kinematical correc- 
tion largely reduces the rate in the high-g^ region. After the 6-quark mass is consistently 
taken into account in the CSS-HQ resummation formalism [63], the position of the peak 
in the qx distribution shifts to a lower value, while the rate in the high-g^ region remains 
unchanged. The combination of the kinematic correction and heavy-quark mass correction 




10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 

q^iggsfCeV) q-,n>9gs(GeV) 



Figure 6: Comparison of q T distributions predicted by TOT(l), PERT(l) and PYTHIA, for Higgs 
boson produced via bb fusion at the Tevatron Run-2 (left) and LHC (right) for ttia =100 and 300 
GeV respectively. 
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makes the resummation predictions resemble more the PYTHIA predictions, as shown in 
Fig. § 

As noted in Sec. ||, we assumed the non-perturbative functions in the CSS-HQ resum- 
mation formalism for the heavy bottom quark are the same as those for the light quarks in 
our numerical calculations. However, additional non-perturbative dynamics, such as that 
associated with "intrinsic" heavy quarks [73], might be present in the heavy- flavor channels. 
In the future, when experimental data (likely, from the associated production of Z boson 
with bottom (anti)-quark [74-76], and from the t-channel single-top production [77-86]) 
becomes available to measure the bottom quark parton distribution function together with 
those non-perturbative functions, we will be able to improve our theoretical prediction on 
the qr distribution of any hard-scattering process that is initiated by bottom quark in- 
teraction in the initial state. Nevertheless, we expect the qualitative comparison between 
TOT(l) and PYTHIA predictions shown in Fig. |6| to stay valid, as these non-perturbative 
variations at impact parameters b > bo/irif, 0.25 GeV -1 would result in mild modula- 
tions in qx space at qx < 30 GeV. More importantly, differences in the shape of the two 
predictions in the large-(?T region have implications for the discovery potential of the Higgs 
boson, as they will affect the significance of the signal event after imposing the kinematic 
cut on the transverse momentum of the Higgs boson to suppress background or to enable 
the reconstruction of the signal kinematics. 




Figure 7: Comparison of the integrated rates, deduced form Fig. 6, as a function of the minimal 
qr value taken in the integration at the Tevatron Run-2 (left) and LHC (right) for tua =100 and 
300 GeV, respectively. 
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